import numpy as np
MOD = 10 ** 9 + 7

def power(A, n):
	res = np.identity(n=A.shape[0])
	for _ in range(n):
		res @= A
	return res

def power2(A, n):
	res = np.identity(n=A.shape[0])
	pow = A.copy()
	while n:
		if n & 1:
			res @= pow
		pow @= pow
		n >>= 1
	return res

def power3(A, n):
	eigenvalues, X = np.linalg.eig(A)
	Λ = np.diag(eigenvalues)
	X_inv = np.linalg.inv(X)
	res = X @ (Λ ** n) @ X_inv
	return res

def power4(A, n):
	return np.linalg.matrix_power(A, n)


if __name__ == "__main__":
	A = np.array([[1,1]
	              [1,0]])
	n = 10
	print(power(A, n))
	print(power2(A, n))
	print(power3(A, n))
	print(power4(A, n))

	for _ in range(100):
		A = np.random.randint(1, 10, size=(3,3))
		n = np.random.randint(1, 10)
		r1 = power(A, n)
		r2 = power2(A, n)
		r3 = power3(A, n)
		r4 = power4(A, n)
		if not np.allclose(r1, r2, r3, r4):
			print('error!')
			print(A)
			print(n)
			print(r1)
			print(r2)
			print(r3)
			print(r4)
			break
	print('done!')